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In this article, Euler-Lagrangian dynamics explain that the two particle interaction has non- 
conservative forces about the frame of the center of mass. This interpretation clarifies the underlying 
interaction and the system descriptions become advantages for MD simulations. 



PACS numbers: 

A Molecular Dynamics (MD) Simulation is based on 
Newton's second law on time evolution within the frame- 
work of classical mechanics. The motion of granular dy- 
namics can be described by Euler-Lagrangian dynam- 
ics. From a comprehensive point of view, the Euler- 
Lagrangian dynamics is presented to clearly and easily 
identify the pair interaction and contacting time among 
the individual particles. The kinetic energy and the po- 
tential energy are used to obtain the Lagrangian of the 
particle with the generalized coordinate and its associ- 
ated granular interaction term. Li this work, the system 
of N particles was the only constraint used to move the 
particles through the external gravitational force and the 
energy dissipated in friction. 

For a granular system of N particles moving in three di- 
mensional space with constraints, D'Alembert's principle 
requires knowledge of the constraint forces, (typically and 
static) , kinetic frictional force and particle sliding. What 
was needed was a description of the system that makes 
compromises when dealing with the constraint relations. 
D'Alembert's principle provides such a description for 
systems involving holonomic or nonholonomic constraints 
in allowing virtual displacements, ^r^ = X]a"'^9a i^^ con- 

a 

trast with a real displacement rfr^ . There are two types of 
constraints. The first one is holonomic constraints which 
can be solved for using kinematics; here there is no work 
(F-'^' • 5vi = 0) because virtual displacements Svi are 

(c) 

orthogonal to the corresponding constraint forces F) ' . 
The second is nonholonomic constraints which must be 
solved for using dynamics; here consequently we have a 
work term (F-^^ x 5vi ^ 0). 

As a consequence of these constraint relations, this 
problem is solved by choosing Lagrange's equation. The 
Euler-Lagrange equations, L{qa,qa) in the generalized 
coordinates, Qa for a discrete system is 



d dL dL 
dt dqa dqa 



= Qa 



(1) 



Here L = T — V — U is the Lagrangian of the system. 

The generalized constraint forces, Qa are derived from 
D'Alembert's principle for the virtual work for applied 
forces. Two nonconservative components are associated 
with respect to qa and qa- The generalized constraint 
forces on the right-hand side are given by 
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dvi 
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(2) 



where F^^ ' is only for the friction force and F is known as 
Rayleigh's dissipation function F = ^X^'^u^a^ (if '^ij — 

a.h 

Cji , the symmetric tensor). 

The Lagrange equation with Rayleigh's dissipation 
function and contact friction becomes 



( f) 
where Qa 



d dL dL 
dt dqa dqa 



dqa 



= Qi'' 



(3) 
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j ■ ■ -Q-^- The generalized constraint force 
is presented by a frictional force. 

In the case in which particles were assumed to be the 
spherical rigid bodies located in the gravitational field, 
the generalized coordinates were decoupled into the cen- 
ter of mass coordinates (the axes of the incrtial coordi- 
nate system) and the Euler angles {9,(j),ip). Since the 
spherical rigid body is symmetrical for the Euler angles, 
the moment of inertia for an ideal homogeneous sphere of 
radius R is I = ^mR^ is about the principal axis of the 
center of mass. From the relations above we can compute 
one term relating the kinetic energy T of the translational 
motion and the spin rotation itself about the center of 
mass. Another term relating the potential energy of the 
gravity potential V and a pairwise interaction potential 
U (rij), is given where the distance between the centers 
of the two particles is r^ = jr^ — rj| with respect to a 
fixed frame and r^, Vj is the relative position vector of 
the center of mass rui , mj . 

The distance between particles i and j about a fixed 
point O, the relative position about the normal and tan- 
gential direction, is independent. Above all, the relative 




where f^ — r 



'^jl^ij S-Hd Tij 



Yi — Tj , the relative velocity 








FIG. 1: Schematic of the vector analysis contacting on particle 
acting j — >■ i 



angular velocity in the rotating frame is contributed to 
the angular displacement fr-om the equation 



Ui X 



Ri 



Hi + R 



-Tij + ojj X ■' r,j = uj,j X r„- (4) 



The relationship of the overlap displacement is defined 



as 



5rij= {vij ■ 5ij) fjj +fy X {6ij X f^) - / usij x Vijdr (5) 



Here the center of angular velocity is uJij = 
[RiiOi + RjOJj) / (Ri + Rj) , with angular velocities uJi and 
ujj in a rotating coordinate frame with each particle at 
the center of mass of the system during contact time tc- 
The relative normal displacement is 

frijil = (f.y ■Sij)rij (6) 

and the relative tangential displacement is 

6nj_i_ = r^j X (% X Tij) - / UJij X fydr (7) 

A comparison with the above gives the relative velocity 



as 



dr. 



dt 



^ ij tj ij 



Vo has the com- 



in term of j . 

The relative velocity f y- = v^ == 
ponents of the relative normal displacement 

VijII = (fjj •V,y)f,y (8) 

and the relative tangential velocity is given by 



fiji^ 



r,-,- X (v,-,- X r. 



i-ij 



'ij 
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(9) 



(fij-Vij)f 



I] J •■I] 



where f^ x (v^j x f jj) = v, 

When granular particles have pairwise collisions, they 
collide instantly and have an infinitesimally small time 
for the Hertz's potential. This impulse imparts inter- 
actions between each particle. Figure shows analytic ge- 
ometry, before instantaneous deformation, at which point 
the bodies touch at an adjacent point. The curvatures 
of the two bodies are represented with radii of curvature 
Ri and Rj. At the contact area, the bodies are com- 
pressed with smooth curvature. The shape of the particle 
is changed by the curvature due to the elastic restoration 
force. Previously, Hertz's theory had been used with a 
quadratic equation. However, contact pressure provides 
displacement under the surface. This contact pressure 
is distributed throughout the curvature. The pressure 
distribution function acts at each contact radius a. The 
distribution of normal stress in the contact area as a func- 
tion of distance was also reported [2] . The granular par- 
ticles behave like hard sphere which contact one another 
with a deformable contact force. The forces between the 
two particles can be approximated well by Hertz theory of 
elasticity. The granular interaction potential is a power 
law from Hertz theory. Hertz contact is presented to 
determine each granular surface interaction for granular 
flows. The generalized expression for the normal contact 
force acting on a granular body, in which the governing 
equation is either loaded in tension or in compression, is 
obtained from the stress and strain relation of the mate- 
rial. The Hertz potential is derived from the results of 
stress-strain analysis that are attributed to granular flow 
sensitivity. Granular particle interactions have a virtual 
depth in the relation with a remarkable normal contact 
of relative deformation in terms of a tangential contact 
by a friction on a contact surface. A substantial solution 
to the problem of a normal contact has been known as 
Hertz's theory. 

We are led to the relations of the overlap displace- 
ment as deflned by 5rij where the normal displacement 
(5rjj|| = Ri + Rj — |ri — r^l is the depth of indentation or 
normal overlap for the contact between the two bodies 
of radii Ri and i?2- When Svij^^ > 0, a contact force is 
generated that otherwise would have zero potential. The 
noncohesive granular interaction is described by Hertz- 
Mindlin potential between particles. 
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(10) 
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FIG. 2: Schematic of the interaction between two particles at contact with position vectors r^, r^: The left figure of a set of a 
overlap shape and the right figure of the deformable shape of Hertz-Mindlin force 



where fc„ = ^E* \fE* is the stiffness of the pair- 
wise interaction proposed in Hertz theory (Hertz 1895) 
and the tangential term kt = %G*y/~R* is the tan- 
gential force (Mindlin 1949) with the effective Ra- 
dius R* — RiRj/{Ri + Rj) and the effective Young's 

modulus E* = fi^ + ^] , E the Young's 
modulus and the effective Shear modulus is G* — 

_\ ^l\ -r ^/ j^ _\ 12\ u \ ^ Poisson s ratio is v for 



Ei ' Ej 

the contacting particle i and j respectively [HIllll]. 

Comparisons of this with the Lagrangian obtained pre- 
viously can be expressed by the gravitational potential 



to the other particle and orthogonal to the constraint 
force. If both particles are at a rest state (/is < /ifc), 
the sphere at contact experiences the only force on the 
pebble from the static friction (/x^, the static coefficient 
of friction) which is perpendicular to the curved part of 
particle. The normal and tangential direction is inde- 
pendent, and then each term must be separate. The 
motion of the ith particle is given by a contact force 
and a damping force during collisions. If the tangen- 
tial force is generated, the generalized constraint force is 



q\I' 



-sign {5rij±) min (fi F^ 



-I) 



ij±\) '■ij±- 



Thus, the Lagrange's equation of the net binary system 
can be expressed in term of Vij 



V (r,) = m,ri ■ g 



(11) 



where g is the gravitational acceleration vector. 

The pair dissipative function can be derived from 
Rayleigh's dissipation function 



jr = -jm*Srl 



(12) 



where 7 is a dissipative constant or viscous damping co- 
efficient given a coefRcient of normal restitution e for a 
contacting time with normal and tangential contact. 

When the particle i is moving to another particle j at 
contacting time, it follows that the force of constraint 
exerted by kinetic friction {fj,k, the kinetic coefRcient of 
friction) of each particle has done real work after the 
collision. Thus the virtual displacement, Sr, is tangent 



L^ -TO* \6hjf -Uirij) 



(13) 



where m* = mijrij/ (rrii + rrij) is the reduced mass of the 
system. 

As a consequence of the decoupling, the Lagrange 
equations of motion separate into two equations, one for 
the center of mass coordinates and another for the Eulcr 
angles. Each set can be analyzed independently to ap- 
ply this prescription to a two body problem for a virtual 
displacement. The present simulations use a normal and 
tangential contact model comprised of the following: 



\Shjf^\6hnf+(^ + ^)\Srr,^\' (14) 



where C — i —^ H V 

\mili rrijlj 



from the particle moment of 



inertia li = li/niiRf and Ij — Ij /rrijR'^ about the torque 
acting on two particle. 

Finally, we can implement Euler-Lagrangian dynam- 
ics to find out the governing equation and contact time 
used to simulate the MD method. Although the theo- 
retical approach is Hertz's contact to the Lagrangian so- 
lution, problems arise including general friction through 
contacting area where Fji — —Fji, Newton 3rd law and 
the tangential component Fij± and 5rijj_ which can be 
described to meet the requirement of Coulomb yield crite- 
rion if |Fijj^| < HgF^ju . Until the particles are separated, 
the potential calculation is performed for the interaction 
only once on each pair of particles. The equations of 
motion follow directly from the Lagrangian formulation 
described above. 

The normal component of the contact force can be 
written as 



Fy||=fc„(5r^^ 



3/2 



7„TO V, 



v\\ 



(15) 



The shear component of the contact force can be writ- 
ten as 



ij_L 



rjji -Jtm Vij_L 



(16) 



The equation of motion will be expressed in normal 



and tangential directions on contact condition, |F 



ij_L 



> 



f^sFijW, Fy_L == -/ifcFij||Vjj_L and jtm*Vij± which gen- 
erates the rolling friction from the center of particle j 
to the center of particle i. This component in the direc- 
tion of motion is attributed by performance of the Hertz- 
Mindlin theory to improve the terms of the tangential 
force. The simulation will demonstrate the essential re- 
sults of the contact deformable interaction. Two body 
forces are represented by the potential, which is enough 
to permit the surface effects. The governing equation, 
constructed from two-body functions, is applied to the 
many body system problem with a contact dynamic us- 
ing a MD simulation (or DEM). 
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